clear all;clc;close all;
load('C:\CloudStation\ResearchNote\DiffusionSpin_NF\PumpTrans\PumpTrans_sendtoDW.mat')
Pin = cell2mat(inputPower);Pin=0.7*Pin;
Pout = cell2mat(outputPower);
T = regexp(TempLegendcell,'\d+(\.)?(\d+)?','match');
T = cell2mat(arrayfun(@(x) str2double(T{x}),1:7,'Uniformoutpu',false));

% density=[2.50,6.01,9.00,12.0,15.0,17.9,21.1]*10^12; %1/cm^3
% sigma=0.5*10^-12; %cm^2
% L=1.5;
% OD = density*sigma*L;
% figure()
% plot(T,[0.92,3.5075,5.55,8.0,10.55,12.625,15.825],'-');hold on;
% plot(T,OD,'-*');
% xlabel('T');
% ylabel('OD')
% legend('Exp','Fit')

%
% Pin=Pin/0.7;
% figure()
% hold on;
% arrayfun(@(x) plot(Pin(:,x),Pout(:,x)./Pin(:,x),'.-'),1:length(T))
% legend(TempLegendcell)
% xlabel('P_{in}');ylabel('P_{out}/P_{in}')
% 
% Prel=1.3;
% abs = 1-Prel*lambertw((Pin/Prel).*exp((Pin/Prel)-OD))./Pin;
% figure()
% hold on;
% arrayfun(@(x) plot(Pin(:,x),1-Pout(:,x)./Pin(:,x),'.'),1:length(T))
% arrayfun(@(x) plot(Pin(:,x),abs(:,x),'-'),1:length(T))
% legend(TempLegendcell)
% xlabel('P_{in}');ylabel('abs')
% 
% figure()
% hold on;
% arrayfun(@(x) plot(Pin(:,x),1-Pout(:,x)./Pin(:,x),'.'),1:length(T))
% arrayfun(@(x) plot(Pin(:,x),abs(:,x),'-'),1:length(T))
% legend(TempLegendcell)
% set(gca,'YScale', 'log','XScale','log');
% xlabel('P_{in}');ylabel('abs')
% %%
OD = [0.92,3.5075,5.55,8.0,10.55,12.625,15.825];
absorption = ((Pin - Pout)./Pin);
rescale_Pin = Pin./OD;

Prel=1.48;
abs_lam = 1-0.975*Prel*lambertw((Pin/Prel).*exp((Pin/Prel)-OD))./Pin;

figure();
hold on;
arrayfun(@(x) plot(rescale_Pin(:,x),1-Pout(:,x)./Pin(:,x),'.-'),3:length(T))
set(gca,'YScale', 'log','XScale','log');
arrayfun(@(x) plot(rescale_Pin(:,x),abs_lam(:,x),'-'),3:length(T))
xlabel('P_{in}/OD');ylabel('abs')
legend(TempLegendcell{3:length(T)})

w=0.24;
abs_int_lam=zeros(length(Pin),length(OD));
for i=1:length(Pin)
    for j=1:length(OD)
    abs_int_lam(i,j) = 1.0 - 0.975*Prel*int_lamber(Pin(i,j)/Prel, OD(j), w)/Pin(i,j);
    end
end

figure();
hold on;
arrayfun(@(x) plot(rescale_Pin(:,x),absorption(:,x),'.-'),3:length(T))
set(gca,'YScale', 'log','XScale','log');
arrayfun(@(x) plot(rescale_Pin(:,x),abs_int_lam(:,x),'-'),3:length(T))
xlabel('P_{in}/OD');ylabel('abs')
legend(TempLegendcell{3:length(T)})

figure();
hold on;
arrayfun(@(x) plot(Pin(:,x),absorption(:,x),'.-'),3:length(T))
% set(gca,'YScale', 'log','XScale','log');
arrayfun(@(x) plot(Pin(:,x),abs_int_lam(:,x),'-'),3:length(T))
xlabel('P_{in}/OD');ylabel('abs')
legend(TempLegendcell{3:length(T)})

function res = int_lamber(Pin, OD, w)
    beamshape = @(x,y) exp(-((x-0.5).*(x-0.5) + (y-0.5).*(y-0.5))/2/w^2)/(2*pi*w^2*erf(1/(2*sqrt(2)*w))^2);
    I_xyz = @(x,y) lambertw(Pin.*beamshape(x,y).*exp(Pin.*beamshape(x,y)-OD));
    res= integral2(I_xyz,0,1,0,1);
end